#!/usr/bin/python2.4

# C6304: copyright missing
# pylint: disable-msg=C6304

# C6409: lexical structure of constants
# pylint: disable-msg=C6409

# {([ -- do not delete this line; it is used for file-spanning brace matching

"""Plans route.

Nots:
  To Capture images under X11, try using:
    xwd | xwdtopnm | tee oor.ppm | ppmquant 256 | ppmtogif > oor.gif

"""

import csv
import math
import optparse
import os
import sys
import time

import dfoimagereader
import j42polars
import matrix

from matrix import DM
from matrix import GreatCircleDistance
from matrix import LatLon
from matrix import LatLonDeg
from matrix import Polar
import timeformat
import Tkinter
import uwimagereader
import vectorfield
import windbarb
import xtide_jdf

# TODO(rrh): Add InfTime as a flag
InfTime = 2000000000


def DoTrace():
  return False


class Conditions(object):
  """Holds observed conditions at a certain time.

  There's an array of these indexed by gridded LatLon.

  These conditions are written to a .csv file by the (slow) planning phase,
  and read during the display phase.
  """

  def __init__(self, model_time=InfTime):
    self._model_time = model_time  # unix epoch time (secs)
    self._come_from = None  # where we sailed from
    self._true_wind = Polar(0, 0)  # true wind relative to bottom
    self._true_tide = Polar(0, 0)  # true tidal current relative to bottom
    self._boat_cow = 0  # course over water (compensates for tide)
    self._boat_cog = 0  # course over ground
    self._boat_sog = 0  # speed over ground
    self._boat_dog = 0  # distance over ground
    self._distinguished = False  # this is a point on the route start..end

  def GetTime(self):
    return self._model_time

  def SetTime(self, model_time):
    self._model_time = model_time
    return self._model_time

  # TODO(rrh): Add more setters/getters and/or use properties.

  def __str__(self):
    return ('time=%10d come_from=%s TW=%s tide=%s '
            'COW=%3.0f COG=%3.0f SOG=%4.2f DOG=%4.2f') % (
                int(round(self._model_time)),
                self._come_from,
                self._true_wind,
                self._true_tide,
                self._boat_cow,
                self._boat_cog,
                self._boat_sog,
                self._boat_dog)

  def CSVFormat(self):
    """Returns the value of this condition as CSV."""
    come_from_ll = self._come_from
    if not come_from_ll:
      inf = float('inf')
      come_from_ll = matrix.LatLon(matrix.DM(inf, inf), matrix.DM(inf, inf))
    fmt_vals = [
        ('%10d', int(round(self._model_time))),
        ('%s', come_from_ll.CSVFormat()),
        ('%s', self._true_wind.CSVFormat()),
        ('%s', self._true_tide.CSVFormat()),
        ('%4.0f', self._boat_cow),
        ('%4.0f', self._boat_cog),
        ('%5.3f', self._boat_sog),
        ('%5.3f', self._boat_dog),
        ('%s', self._distinguished),
        ]
    return ','.join([fmt % val for fmt, val in fmt_vals])

  @staticmethod
  def BuildFromCSV(row, i):
    """Convert fields in a scanned CSV row into an instance."""
    condition = Conditions()
    condition._model_time = float(row[i+0])

    come_from_ll, i = LatLon.BuildFromCSV(row, i+1)
    if come_from_ll.lat > 10000:  # Bogus value, hopefully infinity
      condition._come_from = None
    else:
      condition._come_from = come_from_ll

    condition._true_wind, i = Polar.BuildFromCSV(row, i+1)
    condition._true_tide, i = Polar.BuildFromCSV(row, i+1)

    condition._boat_cow = float(row[i+1])
    condition._boat_cog = float(row[i+2])
    condition._boat_sog = float(row[i+3])
    condition._boat_dog = float(row[i+4])
    condition._distinguished = bool(row[i+5])

    return condition, i+6-1  # the caller assumes the whole item is length 1


class GridDisplayer(object):
  """Drives Tkinter to display the grid of interest.

  Handles button press events, updates display, etc.
  """

  def __init__(self, start_time, boatconditions,
               windvectorfield, tidalvectorfield,
               island, isnearland,
               display_images=True):
    self._display_root = Tkinter.Tk()
    self._start_time = start_time
    self._boatconditions = boatconditions
    self._windvectorfield = windvectorfield
    self._tidalvectorfield = tidalvectorfield
    self._island = island
    self._isnearland = isnearland
    self._display_images = display_images

    self._minlat, self._maxlat, self._minlon, self._maxlon = GridBoundsOrdered()

    self._pixels_per_minute_lat = 6
    self._pixels_per_minute_lon = int(
        round(self._pixels_per_minute_lat * math.cos(math.radians(48.0))))

    self._dlat = self._maxlat - self._minlat
    self._dlon = self._maxlon - self._minlon

    self._height_pixels = int(self._dlat) * self._pixels_per_minute_lat
    self._width_pixels = int(self._dlon) * self._pixels_per_minute_lon

    self._route_display_canvas = Tkinter.Canvas(
        width=self._width_pixels,
        height=self._height_pixels,
        background='white')
    self._route_display_canvas.grid(row=0, column=0)
    self._route_display_canvas.bind('<ButtonPress-1>', self.ButtonPressCallback)

    self._quitbutton = Tkinter.Button(
        width=10,
        height=2,
        text='Quit',
        background='red')
    self._quitbutton.grid(row=1, column=0)
    self._quitbutton.bind('<ButtonPress-1>', self.QuitCallback)

    self._frame0 = Tkinter.Toplevel()
    self._frame0.title('Tide Image from DFO')

    self._tide_image_label = Tkinter.Label(self._frame0)
    self._tide_image_label.grid(row=0, column=1)

    self._frame1 = Tkinter.Toplevel()
    self._frame1.title('Wind Image from UW')
    self._wind_image_label = Tkinter.Label(self._frame1)
    self._wind_image_label.grid(row=0, column=2)

    self.DrawCoastLine(WACoastLine(), 'WA')
    self.DrawCoastLine(BCCoastLine(), 'BC')

    isochron_tagname = 'isochron'
    isochron_color = 'brown'
    isochronpoints = ExtractIsoChrons(self._boatconditions, self._start_time)
    for i in sorted(isochronpoints):
      ll_points = isochronpoints[i]  # TODO(rrh): unsorted
      pixel_points = []
      for ll in ll_points:
        pixel_points.extend(self.LatLonToPixels(ll))
      if len(pixel_points) > 2:
        kwargs = {
            'smooth': 'true',
            'width': 2,
            'fill': isochron_color,
            'tags': [isochron_tagname],
            }
        self._route_display_canvas.create_line(*pixel_points, **kwargs)

  def DrawCoastLine(self, ll_points, name, color='black'):
    """Draw a coastline.

    Args:
      ll_points: A list of LatLon objects making up the coastline.
      name: The name of the coastline.
      color: The color of the coastline.
    """
    pixel_points = []
    for ll in ll_points:
      pixel_points.extend(self.LatLonToPixels(ll))
    kwargs = {
        'smooth': 'false',
        'width': 2,
        'fill': color,
        'tags': [name],
        }
    self._route_display_canvas.create_line(*pixel_points, **kwargs)

  def DrawIsoChrons(self, interval=5, color='blue'):
    """Draw the isochrons.

    Args:
      interval: The interval, in hours, to use between isochrons.
      color: The color to draw the isochrons.
    """
    pass

  def DrawRoute(self, waypoints, color='green'):
    """Draw a route, and labels the time at each waypoint.

    Args:
      waypoints: A list of LatLon objects for the waypoints of the route.
      color: The color of the coastline.
    """
    trace = True

    self._route_display_canvas.delete('route')  # Delete the display of the old route

    pixel_points = []
    for wpt in waypoints:
      pixel_points.extend(self.LatLonToPixels(wpt))
    kwargs = {
        'smooth': 'false',
        'width': 2,
        'fill': color,
        'tags': ['route'],
        }
    self._route_display_canvas.create_line(*pixel_points, **kwargs)

    # Draw the tidal vector field at time of the selected end waypoint.
    end_wpt = waypoints[-1]
    condition = self._boatconditions[end_wpt]
    model_time = condition.GetTime()
    minlat, maxlat, minlon, maxlon = GridBoundsOrdered()

    # Get and update the image for the end waypoint at end_wpt.
    if self._display_images:
      p, tidal_attribution = self._tidalvectorfield.FullIndex(
          model_time, end_wpt)
      p, wind_attribution = self._windvectorfield.FullIndex(
          model_time, end_wpt)

      if tidal_attribution:
        tidal_attribution = tidal_attribution.strip()
        tidal_attribution = tidal_attribution[1:-1]  # strip off quotes
        print 'tidal_attribution %s' % (tidal_attribution,)
        self._tide_image_image = Tkinter.PhotoImage(file=tidal_attribution)
        self._tide_image_label.config(image=self._tide_image_image)

      if wind_attribution:
        wind_attribution = wind_attribution.strip()
        wind_attribution = wind_attribution[1:-1]  # strip off quotes
        print 'wind_attribution %s' % (wind_attribution,)
        self._wind_image_image = Tkinter.PhotoImage(file=wind_attribution)
        self._wind_image_label.config(image=self._wind_image_image)

    # Display wind and tidal fields on 6minute latlon grid;
    # Note that we have data for 1minute latlon grids, but that's too
    # much to display.
    for lat in reversed(xrange(minlat, maxlat + 1, 6)):
      for lon in xrange(minlon, maxlon + 1, 6):
        latlon = matrix.LatLon(lat, lon)
        if InObstacle(self._island, self._isnearland, latlon):
          continue

        s = matrix.TransformStack()
        tx, ty = self.LatLonToPixels(latlon)
        s.Translate(matrix.Point(tx, ty))

        s.Push()  # {
        true_wind = self._windvectorfield.Index(model_time, latlon)
        s.Scale(20.0, -20.0)  # compensate for mirrored X11 Y coordinate system
        VectorWindBarbGenerator(self, s, true_wind, 'orange', rotation=-60.0)
        s.Pop()  # }

        s.Push()  # {
        true_tide = self._tidalvectorfield.Index(model_time, latlon)
        VectorTidalVectorGenerator(self, s, true_tide, 'purple')
        s.Pop()  # }

    # Label the arrival time along some of the waypoints.
    # To reduce clutter in the diagram, only print each wpt_period wpts
    wpt_period = 3
    if wpt_period <= 0:
      wpt_period = 1

    for wpt, wpt_number in zip(waypoints, list(xrange(0, len(waypoints)))):
      if wpt_number % wpt_period == 0:
        condition = self._boatconditions[wpt]
        model_time = condition.GetTime()
        p, tidal_attribution = self._tidalvectorfield.FullIndex(model_time, wpt,
            trace=False)
        p, wind_attribution = self._windvectorfield.FullIndex(model_time, wpt,
            trace=False)

        tx, ty = self.LatLonToPixels(wpt)
        formatted_time = timeformat.FormatTime(model_time)

        if trace:
          print '%2d %s ==> @%s ==> {%s} %s %s' % (
              wpt_number,
              wpt,
              formatted_time,
              condition,
              tidal_attribution,
              wind_attribution)

        # displayed_text = formatted_time + ' ' + str(wpt_number)
        displayed_text = 'w%02d @%s' % (wpt_number, formatted_time[11:])
        self._route_display_canvas.create_text(
            tx, ty,
            anchor='w',
            fill='black',
            text=displayed_text,
            tags=['route'])

        s = matrix.TransformStack()
        s.Translate(matrix.Point(tx, ty))
        s.Push()  # {
        s.Scale(20.0, -20.0)  # compensate for mirrored X11 Y coordinate system
        VectorWindBarbGenerator(self, s, condition._true_wind, 'red', rotation=-60.0)
        s.Pop()  # }

        VectorTidalVectorGenerator(self, s, condition._true_tide, 'blue')

    if trace:
      print ''

  def QuitCallback(self, event):
    """Called on a Tkinter button press event in the quit button."""
    print 'QUIT!'
    sys.exit(1)

  def ButtonPressCallback(self, event):
    """Called on a Tkinter button press event in the canvas.

    Will draw the best route to the point that was clicked.

    Args:
      event: event on the canvas; contains x and y in screen coordinates.
    """
    ll = self.PixelsToLatLon(event.x, event.y)
    if DoTrace():
      print 'ButtonPress: %d %d ==> %s at %g' % (
          event.x,
          event.y,
          ll.FullPrint(),
          self._boatconditions[ll].GetTime()
          )
    waypoints = GoalExtractList(self._boatconditions, ll)
    self.DrawRoute(waypoints)

  def Dump(self):
    print 'bounding lat: (%g, %g), bounding lon: (%g, %g)' % (
        self._minlat, self._maxlat,
        self._minlon, self._maxlon,)
    print 'dlat=%g (minutes)' % (self._dlat,)
    print 'dlon=%g (minutes)' % (self._dlon,)
    print 'pixels_per_minute_lat=%g (pixels)' % (self._pixels_per_minute_lat,)
    print 'pixels_per_minute_lon=%g (pixels)' % (self._pixels_per_minute_lon,)
    print 'height_pixels=%g (pixels)' % (self._height_pixels,)
    print 'width_pixels=%g (pixels)' % (self._width_pixels,)

  def PixelsToLatLon(self, x, y):
    """Given Tkinter/X11 screen coordinates, convert to a LatLon.

    Args:
      x: The Tkinter/X11 x screen position.
      y: The Tkinter/X11 y screen position.
    Returns:
      The corresponding LatLon.
    """
    lat = self._maxlat - y / self._pixels_per_minute_lat
    lon = self._minlon + x / self._pixels_per_minute_lon
    return LatLon(lat, lon)

  def LatLonToPixels(self, ll):
    """Given a LatLon, convert to Tkinter/X11 (x, y) screen coordinates.

    Args:
      ll: A LatLon point.
    Returns:
      An (x, y) pair for the Tkinter/X11 screen coordinates.
    """

    x = ll.lon - self._minlon
    y = -(ll.lat - self._maxlat)  # compensate for X11 coordinates
    x *= self._pixels_per_minute_lon
    y *= self._pixels_per_minute_lat
    if DoTrace():
      print 'Convert %s ==> %g %g' % (ll, x, y,)
    return (x, y)

  def Run(self):
    """Runs the display in the Tk loop."""
    self._display_root.mainloop()


def GridBoundsOrdered():
  """Returns the bounding box.

  Returns:
    A quad (minlat, maxlat, minlon, maxlon).
  """
  bounds = GridBounds()
  minlat = min(bounds[0].lat, bounds[1].lat)
  maxlat = max(bounds[0].lat, bounds[1].lat)
  minlon = min(bounds[0].lon, bounds[1].lon)
  maxlon = max(bounds[0].lon, bounds[1].lon)
  return minlat, maxlat, minlon, maxlon


def MakeBoatConditions():
  """Build the BoatConditions matrix, with one Conditions value in each grid.

  Returns:
    The boat conditions matrix: an array from gridded LatLon to Conditions
  """
  boatconditions = {}  # Maps LatLon to the Conditions observed on optimal route
  minlat, maxlat, minlon, maxlon = GridBoundsOrdered()
  for lat in xrange(minlat, maxlat + 1):
    for lon in xrange(minlon, maxlon + 1):
      latlon = LatLon(lat, lon)
      boatconditions[latlon] = Conditions()
  return boatconditions


def MakeLandConditions():
  """Build the data structures describing the land.

  Returns:
    A pair: The first item is a set of LatLon believed to be on land.
    The second item is a set of LatLon believed to be near land.
  """
  island = set()  # set of LatLons believed to be on land
  isnearland = set()  # set LatLons believed to be near land
  coastlines = CoastLines()
  minlat, maxlat, minlon, maxlon = GridBoundsOrdered()
  for lat in xrange(minlat, maxlat + 1):
    for lon in xrange(minlon, maxlon + 1):
      latlon = LatLon(lat, lon)
      is_land = False
      for coastline_points in coastlines:
        is_land |= PnPoly(coastline_points, latlon)
      if is_land:
        island.add(latlon)
        isnearland.add(latlon)
        for dlat in xrange(-5, 6):
          for dlon in xrange(-5, 6):
            latlon_near = LatLon(lat + dlat, lon + dlon)
            isnearland.add(latlon_near)

  return island, isnearland


def ExtractIsoChrons(boatconditions, start_time,
                     interval_hours=3.0, interval_max=24.0*3):
  """Compute points along isochron lines.

  MET == Mission Elapsed Time
  All computations in this routine are done in units of hours, not seconds.

  Args:
    boatconditions: Array mapping gridded LatLon to Conditions.
    start_time: Starting time of the race.
    interval_hours: How far apart the isochrons should be, hours of mission
      elapsed time.
    interval_max: Approximate maximum isocron value.

  Returns:
    An array of points, indexed by integer counting intervals, where
      each interval is width interval_hours.
  """
  nsteps = int(round(interval_max / interval_hours))
  points = {}
  for i in xrange(0, nsteps):
    points[i] = []
  minlat, maxlat, minlon, maxlon = GridBoundsOrdered()
  for lat in xrange(minlat, maxlat + 1):
    for lon in xrange(minlon+1, maxlon + 1):
      this_point = LatLon(lat, lon-0)
      last_point = LatLon(lat, lon-1)
      this_met = (boatconditions[this_point]._model_time - start_time) / 3600
      last_met = (boatconditions[last_point]._model_time - start_time) / 3600
      if this_met <= interval_max and last_met <= interval_max:
        for i in xrange(0, nsteps):
          met = i * interval_hours
          interpolate_lon = None
          if last_met <= met and met < this_met:
            interpolate_lon = matrix.Ratio(
                last_met, met, this_met,
                last_point.lon, this_point.lon)
          elif this_met <= met and met < last_met:
            interpolate_lon = matrix.Ratio(
                this_met, met, last_met,
                this_point.lon, last_point.lon)
          if interpolate_lon:
            isochron_point = LatLon(lat, interpolate_lon)
            points[i].append(isochron_point)

  # Use a stupid O(N^2) algorithm to sort the points on an isochron.
  # Choose the point that is nearest to either end of the growing curve,
  # and place that point next to the end of the curve that it is closest to.
  # Doing this sort prevents odd behaviour with curves that wander around
  # our lines of traversal.
  for i in sorted(points):
    all_points = points[i]
    if not all_points:
      continue
    tovisit = set(all_points)
    end0 = all_points[0]
    tovisit.discard(end0)
    ordered_list = [end0]
    while tovisit:
      best_dist = float('inf')
      for point in tovisit:
        dist0, _ = GreatCircleDistance(ordered_list[0], point)
        dist1, _ = GreatCircleDistance(ordered_list[-1], point)
        if dist0 < best_dist:
          best_dist = dist0
          best_point = point
          best_end = 0
        if dist1 < best_dist:
          best_dist = dist1
          best_point = point
          best_end = 1
      tovisit.discard(best_point)
      if best_end == 1:
        ordered_list.append(best_point)
      else:
        ordered_list.insert(0, best_point)
    points[i] = ordered_list
  return points


def PrintConditionsAsCSV(boatconditions, output_file=sys.stdout):
  """Prints the boatconditions values as a csv file to output_file."""
  for latlon in sorted(boatconditions):
    print >>output_file, '%s,' % (latlon.CSVFormat(),),
    condition = boatconditions[latlon]
    print >>output_file, '%s' % (condition.CSVFormat(),),
    print >>output_file, ''


def ReadConditionsAsCSV(input_file_name):
  """Read from a CSV file all the conditions computed previously.

  Args:
    input_file_name: Name of csv file to read with previously saved
      Conditions.

  Returns:
    An array indexed by gridded LatLon to the observed Conditions
      on the optimal route.
  """
  boatconditions = {}  # Maps LatLon to the Conditions observed on optimal route
  input_fd = open(input_file_name, 'rb')
  reader = csv.reader(input_fd)
  for row in reader:
    latlon, i = LatLon.BuildFromCSV(row, 0)
    condition, i = Conditions.BuildFromCSV(row, i + 1)
    boatconditions[latlon] = condition
    # print 'read Conditions of %s %s' % (latlon, condition,)
  input_fd.close()
  return boatconditions


def PrintConditionsAsGrid(boatconditions):
  """Prints the boatconditions values on a grid, in lat/lon coordinates."""
  minlat, maxlat, minlon, maxlon = GridBoundsOrdered()

  print '%6s' % ('',),
  for lon in xrange(minlon, maxlon + 1):
    print '%4s ' % (matrix.LonMinutesAsString(lon),),
  print ''

  for lat in reversed(xrange(minlat, maxlat + 1)):
    print '%6s' % (matrix.LatMinutesAsString(lat),),
    for lon in xrange(minlon, maxlon + 1):
      point = matrix.LatLon(lat, lon)
      if boatconditions[point]._distinguished >= 0:
        print '*%03d ' % (boatconditions[point]._distinguished,),
      else:
        print '%4.1f ' % (boatconditions[point].GetTime(),),
    print ''


def DummyWindVector(start_time, model_time, unused_latlon):
  """Returns true wind speed (kts) and true wind angle (degrees).

  TODO(rrh): Turn this into a mock, as a subclass of VectorField.

  Args:
    start_time: The starting time of the race (secs).
    model_time: The unix epoch time (secs).
    unused_latlon: The position (LatLon).

  Returns:
    The true wind information, as a Polar.
  """

  delta_t = model_time - start_time  # Time since race began (secs).
  delta_t_wind = 36 * 3600  # Time for wind to reach end_speed (secs).

  ratio = (1.0 * delta_t) / delta_t_wind

  start_true_wind_speed = 25.0
  end_true_wind_speed = 5.0

  true_wind_speed = (
      start_true_wind_speed -
      (start_true_wind_speed - end_true_wind_speed) * ratio)

  if true_wind_speed < 0:
    true_wind_speed = 0.0

  start_true_wind_angle = 330.0
  end_true_wind_angle = 180.0
  computed_true_wind_angle = (
      start_true_wind_angle -
      (start_true_wind_angle - end_true_wind_angle) * ratio)
  tw = Polar(true_wind_speed, computed_true_wind_angle)
  return tw


class VectorWindBarbGenerator(windbarb.WindBarbGenerator):

  def DrawLine(self, handle, p0, p1, color='brown'):
    pixel_points = [p0.x, p0.y, p1.x, p1.y]
    kwargs = {
        'smooth': 'false',
        'width': 2,
        'fill': color,
        'tags': ['route'],
        }
    handle._route_display_canvas.create_line(*pixel_points, **kwargs)


class VectorTidalVectorGenerator(windbarb.TidalVectorGenerator):

  def DrawLine(self, handle, p0, p1, color='blue'):
    pixel_points = [p0.x, p0.y, p1.x, p1.y]
    kwargs = {
        'smooth': 'false',
        'width': 2,
        'fill': color,
        'tags': ['route'],
    }
    handle._route_display_canvas.create_line(*pixel_points, **kwargs)


def BoatSpeed(optimal_sailing, unused_p, unused_model_time, bearing, true_wind):
  """Compute the boat speed as a function of true wind speed and true angle.

  The point p and the model_time are passed in, but we would not expect
  the boat speed to vary as a function of location or of time of day,
  although in reality the chop/waves vary diurnally, as well as the
  ability of the sailors to sail the boat optimally.

  TODO(rrh): compensate for chop slowing the boat down.
  TODO(rrh): compensate for set due to current.

  Args:
    optimal_sailing: A class modelling boat's optimal sailing conditions.
    unused_p: The latlon where we are querying for boat speed.
    unused_model_time: The unix epoch time, secs.
    bearing: Direction the boat is sailing, degrees, relative to the water.
    true_wind: A Polar with the true wind speed of the wind, in kts, and
      true wind direction, degrees, relative to the bottom.

  Returns:
    Boat speed, in kts.
  """
  true_wind_speed = true_wind.Magnitude()
  true_wind_angle = true_wind.Angle()

  bow_wind_angle = abs(true_wind_angle - bearing)

  # Assume that port tack and starboard tack behave identically.
  # Assume that the ability to sail optimally is not a function of time of day.
  if bow_wind_angle >= 180:
    bow_wind_angle = 360 - bow_wind_angle

  tw_relative_bow = Polar(true_wind_speed, bow_wind_angle)
  optimal_boat_speed = optimal_sailing.Index(tw_relative_bow)
  return optimal_boat_speed


def SimulatedSailing(tw):
  """A simple model of sailing speed through the water.

  Args:
    tw: A true_wind polar with the
      true wind speed, relative to the water,
      and true wind angle, relative to the bow of the boat.

  Returns:
    The boat speed, in kts.
  """
  true_wind_speed = tw.Magnitude()
  bow_wind_angle = tw.Angle()
  scale = true_wind_speed / 15.0
  if scale > 1.0:
    scale = 1.0

  if bow_wind_angle > 45:
    if bow_wind_angle < 80:
      boat_speed = scale * 6.5  # pointing
    elif bow_wind_angle < 160:
      boat_speed = scale * 7.5  # reach or spinnaker
    else:
      boat_speed = scale * 7.0  # dead down wind
  elif bow_wind_angle > 35:
    # pinching
    boat_speed = scale * 3.0
  else:
    # in irons
    boat_speed = scale * 0.0
  return boat_speed


def InObstacle(island, unused_isnearland, latlon):
  return latlon in island


def NearObstacle(unused_island, isnearland, latlon):
  return latlon in isnearland


def PnPoly(poly, test):
  """See if the test LatLon is inside of the polygon poly.

  From http://www.visibone.com/inpoly/inpoly.c.txt
  Converted by rrh to Python from the sample c code.
  Assumes that there's no distortion due to latitude's effect on longitude.

  Args:
    poly: A list of points making up the edges of the polygon.
      TODO(rrh): Not sure if the poly needs to be closed or not.
    test: The latlon to test if it is inside of the polygon.

  Returns:
    True or False.
  """
  xt = test.lon
  yt = test.lat
  npoints = len(poly)
  inside = False

  if npoints < 3:
    return False
  xold = poly[npoints-1].lon
  yold = poly[npoints-1].lat
  for i in xrange(npoints):
    xnew = poly[i].lon
    ynew = poly[i].lat
    if xnew > xold:
      x1 = xold
      x2 = xnew
      y1 = yold
      y2 = ynew
    else:
      x1 = xnew
      x2 = xold
      y1 = ynew
      y2 = yold
    if (xnew < xt) == (xt <= xold) and (yt-y1)*(x2-x1) < (y2-y1)*(xt-x1):
      inside = not inside
    xold = xnew
    yold = ynew

  return inside


class RoutePlanner(object):
  """Plans routes given various vector fields, boundaries and conditions."""

  def __init__(self,
               windvectorfield, tidevectorfield,
               starting_name, ending_name,
               start_time, optimal_boat_sailing,
               max_iteration_count=1000):
    self._windvectorfield = windvectorfield
    self._tidevectorfield = tidevectorfield
    self._optimal_boat_sailing = optimal_boat_sailing
    self._start_time = start_time
    self._boatconditions = MakeBoatConditions()
    self._island, self._isnearland = MakeLandConditions()
    self._starting_name = starting_name
    self._ending_name = ending_name
    self._starting_point = WellKnownLocation(starting_name)
    self._ending_point = WellKnownLocation(ending_name)
    self._max_iteration_count = max_iteration_count

  def PlanAndSaveRoute(self, output_file):
    """Plan the route and save results to output_file."""
    print 'Sail From: %15s at %s' % (self._starting_name, self._starting_point)
    print 'Sail   To: %15s at %s' % (self._ending_name, self._ending_point)
    self._grid_bounds = GridBounds()

    print 'Sail   In: %s .. %s' % (self._grid_bounds[0], self._grid_bounds[1])

    self._boatconditions[self._starting_point].SetTime(self._start_time)

    self._starting_set = set()
    self._starting_set.add(self._starting_point)

    self._SailGrid(self._starting_set)

    # Mark the Conditions with the optimal waypoints.
    self._GoalExtract(0, self._ending_point)

    # Save the Conditions to a file for future consumption.
    out_fd = open(output_file, 'w')
    PrintConditionsAsCSV(self._boatconditions, out_fd)
    out_fd.close()

  def _SailTo(self, model_time, p1, p2):
    """Sail from p1 to p2, starting at model_time.

    The wind and tidal vectors valid for p1 is used for the entire segment
    from p1 to p2.

    Args:
      model_time: Unix epoch time (secs).
      p1: starting point (LatLon).
      p2: ending point (LatLon).

    Returns:
      True if getting to p2 from p1 beats all other
        candidate routes to get to p2.
      False if this candidate didn't beat others.
    """
    if InObstacle(self._island, self._isnearland, p2):
      return False

    if NearObstacle(self._island, self._isnearland, p2):
      crosses = CheckIfRouteCrossesCoast(p1, p2)
      if crosses:
        return False

    distance, bearing = GreatCircleDistance(p1, p2)

    true_wind = self._windvectorfield.Index(model_time, p1)
    true_tide = self._tidevectorfield.Index(model_time, p1)

    # boat_speed is speed over water.
    boat_speed = BoatSpeed(
        self._optimal_boat_sailing, p1, model_time, bearing, true_wind)

    if boat_speed:
      delta_t_hours = distance / boat_speed
    else:
      delta_t_hours = float('Inf')

    bearing_cow = bearing
    delta_t_hours, bearing_cow = self._IterateTidalCompensation(
        model_time, p1,
        distance, bearing,
        true_wind, true_tide, delta_t_hours)

    delta_t_secs = delta_t_hours * 3600
    arrival_time = model_time + delta_t_secs

    if (p2 in self._boatconditions and
        arrival_time < self._boatconditions[p2].GetTime()):
      conditions = self._boatconditions[p2]
      conditions.SetTime(arrival_time)
      conditions._come_from = p1
      conditions._boat_sog = boat_speed
      conditions._boat_cow = bearing_cow  # course over water
      conditions._boat_cog = bearing
      conditions._boat_dog = distance
      conditions._true_wind = true_wind
      conditions._true_tide = true_tide
      return True

    return False

  def _IterateTidalCompensation(self,
                                model_time, p1,
                                distance, bearing,
                                true_wind, true_tide, delta_t_hours):
    # Iterate several times to compensate for the tidal offset.
    # We're trying to get a better estimate for delta_t_hours
    # by sailing to a modified point route_vector_offset
    # so that powering through the water from p1 to route_vector_offset,
    # plus tidal offset over the ground from route_vector_offset
    # to p2 yields the true passage time over the bottom from p1 to p2.
    # In empirical practice, 3 iterations or fewer gets us to to
    # within 1 part per 1000.
    trace = False
    for i in xrange(0, 5):
      route_vector = Polar(distance, bearing)
      tide_distance_nm = true_tide.Mag() * delta_t_hours
      tide_distance_polar = Polar(tide_distance_nm, true_tide.Angle())
      route_vector_offset = route_vector - tide_distance_polar
      if trace:
        print (
            '%d delta_t_hours=%5.2f route_vector=%s '
            'tide_distance_nm=%4.2f, tide_distance_polar=%s '
            'route_vector_offset=%s') % (
                i,
                delta_t_hours,
                route_vector,
                tide_distance_nm,
                tide_distance_polar,
                route_vector_offset)
      bearing_cow = route_vector_offset.Angle()
      new_boat_speed = BoatSpeed(
          self._optimal_boat_sailing, p1, model_time, bearing_cow, true_wind)
      if new_boat_speed:
        new_delta_t_hours = route_vector_offset.Mag() / new_boat_speed
      else:
        new_delta_t_hours = float('inf')

      time_difference_hours = new_delta_t_hours - delta_t_hours
      if trace:
        print (
            '    %d %7.4f bearing_cow=%4.0f delta_t_hours=%7.4f '
            'new_delta_t_hours=%7.4f') % (
                i,
                time_difference_hours,
                bearing_cow,
                delta_t_hours,
                new_delta_t_hours)
      if new_delta_t_hours == float('inf'):
        break
      if int(1000*abs(new_delta_t_hours - delta_t_hours)) < 1:
        break
      delta_t_hours = new_delta_t_hours

    if trace:
      print ''
    return delta_t_hours, bearing_cow

  def _SailFrom(self, model_time, p1):
    """Sail from p1 to all the ~120 nearby grid points.

    Args:
      model_time: The unix epoch time that we depart p1 (secs).
      p1: The originating point we are sailing from.

    Returns:
      A set of new bested candidates visited from p1.
    """
    work_set = set()
    for dlat in xrange(-5, 6):
      for dlon in xrange(-5, 6):
        p2 = LatLon(p1.lat + dlat, p1.lon + dlon)
        betters_time = self._SailTo(model_time, p1, p2)
        if betters_time:
          work_set.add(p2)
    return work_set

  def _SailGrid(self, starting_set):
    """Run the simulation to sail the grid.

    This is a classical relaxation algorithm, "repeat until no change".

    Args:
      starting_set: a set of LatLon points representing the starting points to
        evaluate.
    """
    iteration_count = 0
    print 'max_iteration_count=%d' % (self._max_iteration_count,)
    while starting_set and iteration_count < self._max_iteration_count:
      new_starting_set = set()
      iteration_count += 1
      print >>sys.stdout, 'Start step %3d: len(starting_set) = %4d' % (
          iteration_count, len(starting_set))
      sys.stdout.flush()
      for starting_point in starting_set:
        starting_time = self._boatconditions[starting_point].GetTime()
        if DoTrace():
          print 'starting_point=%s at %g' % (starting_point, starting_time,)
          sys.stdout.flush()
        new_starting_set |= self._SailFrom(
            self._boatconditions[starting_point].GetTime(),
            starting_point)
      starting_set = new_starting_set

  def _GoalExtract(self, depth, p_end):
    """Mark the boat_conditions array with the integer waypoint number.

    TODO(rrh): This code may be deprecated.

    Args:
      depth: The waypoint number (perhaps reversed).
      p_end: The ending point of the journey.
    """
    if DoTrace():
      print 'GoalExtract %s' % (p_end,)
      sys.stdout.flush()

    if p_end:
      self._boatconditions[p_end]._depth = depth
    if p_end in self._boatconditions:
      newpoint = self._boatconditions[p_end]._come_from
      self._GoalExtract(depth + 1, newpoint)


def GoalExtractList(boatconditions, p_end):
  """Return latlong waypoints along the way to the goal p_end.

  Do this by tracing back from p_end through the grid of optimal sailings.

  Args:
    boatconditions: An array from gridded LatLon to Conditions.
    p_end: The ending (goal) point of the journey, a LatLon.

  Returns:
    The waypoint list; the first item is the start point,
      and the last item is p_end.
  """
  waypoints = []
  GoalExtractListWorker(boatconditions, waypoints, p_end)
  waypoints.reverse()
  return waypoints


def GoalExtractListWorker(boatconditions, accumulating_list, p_end):
  """Accumulate latLong waypoints along the way to the goal p_end.

  This is a recursive procedure, working from p_end to the starting point.
  Every grid cell is marked with a field _come_from, which ultimately take us
  to the starting point of the calculation.

  Args:
    boatconditions: An array from gridded LatLon to Conditions.
    accumulating_list: A list to append LatLon waypoints, in order.
    p_end: The ending (goal) point of the journey, a LatLon.
  """
  if p_end:
    accumulating_list.append(p_end)
  if p_end in boatconditions:
    condition = boatconditions[p_end]
    if DoTrace():
      print 'GoalExtractList %s ==> %s' % (p_end, condition)
      sys.stdout.flush()
    newpoint = condition._come_from
    GoalExtractListWorker(boatconditions, accumulating_list, newpoint)


def GridBounds():
  """Return a list of LatLon points making up the sailing bounding box."""
  grid_bounds = [
      LatLon(
          DM(48, 40),  # well north of Swiftsure Bank
          -DM(125, 10)  #  well west of furthest point (TODO(rrh): ??)
          ),
      LatLon(
          #ColumbiaBar().lat - 6,  # In minutes
          #Victoria().lon + 6  # In minutes west
          DM(46, 15-6),  # 3 miles south of Columbia bar
          -DM(123, 21-6)  # 6 minutes east of victoria
          ),
      ]
  return grid_bounds


def CheckIfRouteCrossesCoast(p1, p2):
  """Check if the route from p1 to p2 crosses any part of the coast."""
  # TODO(rrh): inefficient, since many terms related to p1 and p2 are common.
  for coastline in CoastLines():
    p3 = coastline[0]
    for p4 in coastline[1:]:
      intersect = matrix.LatLonIntersect(p1, p2, p3, p4)
      if intersect:
        if False:
          print 'Intersect: <%s %s> with <%s %s>' % (p1, p2, p3, p4)
        return True
      p3 = p4
  return False


def SwiftsureBank():
  if False:
    # This point is close to the end of the Boundary line drawn by UW WRF
    return LatLon(48.498*60, -124.745*60)  # decimal degrees to minutes
  else:
    return LatLon(DM(48, 27.0), -DM(124, 35.0))  # From xtide JDF Entrance


def GraysHarborEntrance():
  return LatLon(46.905246*60, -124.174347*60)  # From maps.google.com


def RaceRocksPortAngeles():
  return LatLon(48.230162*60, -123.5289*60)  # From maps.google.com


def ColumbiaBar():
  # return LatLon(DM(46, 15), -DM(124, 2))
  return LatLon(DM(46, 12), -DM(124, 8))  # Approx, from maps.google.com
  # return LatLonDeg(46.205973, -124.136238)  # Approx, from maps.google.com


def Flattery():
  return LatLon(DM(48, 23), -DM(124, 44))  # Flattery (or Tatoosh?)


def WayOffColumnbiaBar():
  return LatLon(DM(46, 15), -DM(124, 52))


def Victoria():
  return LatLon(DM(48, 30), -DM(124, 12))


_NAMES_TO_LOCATIONS = {
    'Victoria': Victoria(),
    'Bar': ColumbiaBar(),
    'ColumbiaBar': ColumbiaBar(),
    'Flattery': Flattery(),
    'SwiftSure': SwiftsureBank(),
    'SwiftShure': SwiftsureBank(),
    'GraysHarbor': GraysHarborEntrance(),
    }


def WellKnownLocation(name):
  return _NAMES_TO_LOCATIONS[name]


def CoastLines():
  return [
      WACoastLine(),
      BCCoastLine(),
      ]


def BCCoastLine():
  """Returns the WA Coastline points, including interior points for fat poly."""
  # Cut and Paste from google maps (right click on a point of interest)
  return [
      LatLonDeg(48.723585, -125.106812),
      LatLonDeg(48.540250, -124.503937),
      LatLonDeg(48.574790, -124.416046),
      LatLonDeg(48.532976, -124.459991),
      LatLonDeg(48.376321, -123.918915),
      LatLonDeg(48.378145, -123.863983),
      LatLonDeg(48.356249, -123.821411),
      LatLonDeg(48.354424, -123.748627),
      LatLonDeg(48.357162, -123.717041),
      LatLonDeg(48.337995, -123.587952),
      LatLonDeg(48.309688, -123.559113),
      LatLonDeg(48.383618, -123.475342),
      LatLonDeg(48.433757, -123.445129),
      LatLonDeg(48.417353, -123.401184),
      LatLonDeg(48.415530, -123.303680),
      LatLonDeg(48.452887, -123.263855),
      LatLonDeg(48.781533, -123.101807),
      ]


def WACoastLine():
  """Returns the WA Coastline points, including interior points for fat poly."""
  return [
      LatLonDeg(46.042736, -123.947754),
      LatLonDeg(46.225453, -123.997192),
      LatLonDeg(46.198844, -123.750000),
      LatLonDeg(46.267240, -123.815918),
      LatLonDeg(46.278631, -124.079590),
      LatLonDeg(46.634351, -124.074097),
      LatLonDeg(46.487047, -123.958740),
      LatLonDeg(46.721035, -123.865356),
      LatLonDeg(46.736096, -124.096069),
      LatLonDeg(46.901492, -124.145508),
      LatLonDeg(46.897739, -124.057617),
      LatLonDeg(46.965259, -123.870850),
      LatLonDeg(47.040182, -124.123535),
      LatLonDeg(46.950262, -124.145508),
      LatLonDeg(47.368594, -124.321289),
      LatLonDeg(47.813155, -124.486084),
      LatLonDeg(47.894248, -124.656372),
      LatLonDeg(48.155093, -124.733276),
      LatLonDeg(48.312428, -124.678345),
      LatLonDeg(48.394106, -124.739456),
      LatLonDeg(48.387722, -124.637146),
      LatLonDeg(48.387722, -124.637146),
      LatLonDeg(48.383618, -124.600754),
      LatLonDeg(48.266741, -124.300003),
      LatLonDeg(48.258970, -124.297943),
      LatLonDeg(48.263541, -124.249191),
      LatLonDeg(48.176618, -124.047318),
      LatLonDeg(48.166543, -123.710861),
      LatLonDeg(48.139974, -123.431396),
      LatLonDeg(48.181654, -123.103180),
      LatLonDeg(47.816843, -122.448120),  # Closure traversal
      LatLonDeg(46.157005, -122.332764),
      LatLonDeg(45.989329, -123.925781),
      ]


def main():
  p = optparse.OptionParser()
  p.add_option('--start', type='str', default='2011-05-12 0930',
               help='When the race is to start.')
  p.add_option('--plan', default=False, action='store_true',
               help='Plan the route, saving results to an intermediate file.')
  p.add_option('--mock', default=False, action='store_true',
               help='Use mock vectorfields for wind, tide and boat polars.')
  p.add_option('--display', default=False, action='store_true',
               help='Display the optimal sailing course for selected point.')
  p.add_option('--display_images', default=True, action='store_true',
               help='Display the tidal and wind images if --display is given.')
  p.add_option('--trace_level', type='int', default=0,
               help='controls some aspects of tracing')
  p.add_option('--max_iteration_count', type='int', default=1000,
               help='maximum number of iterations in the relaxation phase.')
  p.add_option('--ss', default=False, action='store_true',
               help='Report on wind at various major landmarks.')
  p.add_option('--comparetidal', default=False, action='store_true',
               help='compare tidal data from DFO with that from xtide.')

  options, unused_args = p.parse_args()

  start_time = ConvertTimePhrase(options.start)
  print 'Starting time is %s' % (timeformat.FormatTime(start_time),)

  conditions_file_name = 'conditions.csv'

  if options.ss:
    # Report on wind at various easy to eyeball locations,
    # to cross correlate what the map says with what was scraped from the map.
    windvectorfield = uwimagereader.UWWRFVectorField('uw_wrf')
    windvectorfield.ReadFrom(uwimagereader.GetSummaryCSVFileName())
    for ll, name in [
        (SwiftsureBank(), 'SS'),
        (GraysHarborEntrance(), 'GH'),
        (RaceRocksPortAngeles(), 'RRPA'),
        ]:
      for dt in xrange(0, 29):
        delta_t = dt * 3 * 3600
        probe_time = start_time + delta_t
        speed, angle = windvectorfield.Index(probe_time, ll)
        print '%s: %s: %d %02d %4.0f %4.0f' % (
            name,
            timeformat.FormatTime(probe_time),
            probe_time,
            dt * 3,
            speed, angle)

  if options.comparetidal:
    # Compare tidal data scraped from DFO with that scraped from xtide.
    start_time = timeformat.ScanTimeString('2011-05-11 1800')  # start of dfo
    start_time = timeformat.ScanTimeString(options.start)

    vf_dfo = dfoimagereader.DFOVectorField('jdf_dfo')
    vf_dfo.ReadFrom(dfoimagereader.GetSummaryCSVFileName())

    # Requires that extract_xtide.py be run first
    vf_xtide = vectorfield.VectorField('jdf_xtide')
    vf_xtide.ReadFrom('xtide_jdf.csv')

    for hour_step in xrange(0, 3*24):
      model_time = start_time + hour_step * 1 * 3600
      latlon = SwiftsureBank()
      dfo_current = vf_dfo.Index(model_time, latlon)
      xtide_current = vf_xtide.Index(model_time, latlon)
      print '%s %s dfo=%s xtide=%s' % (
          timeformat.FormatTime(model_time),
          latlon,
          dfo_current,
          xtide_current)
    sys.exit(1)

  if options.mock:
    windvectorfield = ConstantWindMock()
    tidevectorfield = ConstantTidalMock()
    optimal_boat_sailing = MotorBoatMock()
  else:
    windvectorfield = uwimagereader.UWWRFVectorField('uw-wrf')
    windvectorfield.ReadFrom(uwimagereader.GetSummaryCSVFileName())
    tidevectorfield = dfoimagereader.DFOVectorField('jdf_dfo')
    tidevectorfield.ReadFrom(dfoimagereader.GetSummaryCSVFileName())
    optimal_boat_sailing = j42polars.J42OptimalSailingComputation()

  if options.plan:
    # Plan the route, writing results to the conditions_file_name
    routeplanner = RoutePlanner(
        windvectorfield,
        tidevectorfield,
        'ColumbiaBar', 'Victoria', start_time,
        optimal_boat_sailing,
        max_iteration_count=options.max_iteration_count)
    routeplanner.PlanAndSaveRoute(conditions_file_name)

  if options.display:
    # Display the route, reading results from the conditions_file_name

    island, isnearland = MakeLandConditions()

    # Read the boat conditions csv file, which summarizes the
    # status of the boat at each grid cell for optimal passage
    boatconditions = ReadConditionsAsCSV(conditions_file_name)
    displayer = GridDisplayer(start_time, boatconditions,
                              windvectorfield, tidevectorfield,
                              island, isnearland,
                              options.display_images)
    if DoTrace():
      displayer.Dump()
    displayer.Run()


def ConvertTimePhrase(name):
  if name == 'now':
    # Start within the hour
    name_to_time = (int(time.time()) / 3600) * 3600
  elif name == 'today':
    # Start at 0900 today
    name_to_time = (int(time.time()) / 86400) * 86400 + 9 * 3600
  else:
    name_to_time = timeformat.ScanTimeString(name)
  return name_to_time


class MotorBoatMock(j42polars.OptimalSailingComputation):
  """Manages optimal sailing[sic] performance data for a perfect motorboat.

  This is used as a mock for J42OptimalSailingComputation.
  """

  def __init__(self):
    super(MotorBoatMock, self).__init__()
    self._complained = False
    self._constant_speed = 10.0

  def _LoadPerformanceItems(self):
    pass

  def Index(self, tw, trace=False):
    """Returns the target boat speed given a true wind angle."""
    if not self._complained:
      print 'WATCHOUT: Using MotorBoatMock with %s' % (self._constant_speed)
      self._complained = True
    if trace:
      print 'MotorBoat tw=%s' % (tw,)
    return self._constant_speed


class ConstantTidalMock(object):
  """Mock up a constant tidal vector."""

  def __init__(self):
    super(ConstantTidalMock, self).__init__()
    self._complained = False
    self._constant_tidal = Polar(10.0, 90)  # current flowing R->L (W->E)

  def FullIndex(self, index_time, latlon, trace=False):
    if not self._complained:
      print 'WATCHOUT: Using ConstantTidalMock with %s' % (
          self._constant_tidal,)
      self._complained = True
    if trace:
      print 'ConstantTidalMock time=%d %s at %s' % (
          index_time, timeformat.FormatTime(index_time), latlon)
    return self._constant_tidal, None


class ConstantWindMock(object):
  """Mock up a constant wind vector."""

  def __init__(self):
    super(ConstantWindMock, self).__init__()
    self._complained = False
    self._constant_wind = Polar(10, 0.0)  # Wind out of due north

  def FullIndex(self, index_time, latlon, trace=False):
    if not self._complained:
      print 'WATCHOUT: Using ConstantWindMock with %s' % (self._constant_wind)
      self._complained = True
    if trace:
      print 'ConstantWindMock time=%d %s at %s' % (
          index_time, timeformat.FormatTime(index_time), latlon)
    return self._constant_wind, None


if __name__ == '__main__':
  main()

# })] -- do not delete this line; it is used for file-spanning brace matching
